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Preface 


This  report  documents  results  of  an  evaluation  of  24  state-of-the-art 
numerical  models  suitable  for  floating  breakwater  applications.  The  work  was 
sponsored  through  funds  provided  to  the  Coastal  Engineering  Research  Center 
(CERC),  US  Army  Engineer  Waterways  Experiment  Station  (WES)  by  the  Civil 
Works  Directorate,  Office,  Chief  of  Engineers  (OCE),  US  Army  Corps  of  Engineers, 
under  the  Coastal  Structures  Evaluation  and  Design  Program  Work  Unit  31679 
titled  "Design  of  Floating  Breakwaters."  OCE  point  of  contact  was  Mr.  Jesse 
Pfeiffer,  Jr. 

The  report  was  prepared  under  WES  Purchase  Order  No.  DACW39-85-M-438l^- 
Mr.  Kozo  Bando,  an  employee  of  Kajima  Corporation,  Japan,  performed  this  study 
as  the  senior  author  while  serving  as  a  resident  engineer  at  TEKMARINE,  Inc. 

The  principal  investigator  of  the  evaluation  was  Dr.  Choule  J.  Sonu,  TEKMARINE, 
Inc.  Mr.  Toshihiko  Takahashi  of  Kajima  Corporation,  Japan,  kindly  provided 
useful  information  to  support  the  study.  Dr.  H.  S.  Chen,  Research  Division, 
CERC,  offered  helpful  technical  comments  during  discussions  of  evaluation 
results. 

This  effort  was  coordinated  by  Mr.  Peter  J.  Grace,  Research  Hydraulic 
Engineer,  Wave  Research  Branch,  Wave  Dynamics  Division  ( WDD ) ,  under  general 
supervision  of  Dr.  James  R.  Houston  and  Mr.  Charles  C.  Calhoun,  Jr.,  Chief 
and  Assistant  Chief,  CERC,  respectively,  and  under  direct  supervision  of  Messrs. 
C.  E.  Chatham,  Chief,  WDD,  and  D.  D.  Davidson,  Chief,  Wave  Research  Branch. 

During  preparation  of  this  report,  COL  Allen  F.  Gram,  USA,  was  Director 
of  WES.  During  publication  of  this  report,  COL  Dwayne  G.  Lee,  CE,  was  Commander 
and  Director.  Dr.  Robert  W.  Whalin  was  Technical  Director. 
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NUMERICAL  MQGB£ 

FOR  A  FLOATING  BREAKWATER 


I.  Scope  of  Study 

The  objective  of  this  study  was  to  identify  an  efficient  numerical 
model  which  will  analyze  and  compute  wave  attenuation  and  moorinq  forces 
related  to  a  floating  breakwater.  Of  particular  interest  was  the  deter¬ 
mination  of  wave  height  attenuation  in  the  lee  of  a  floating  breakwater 
in  a  random  sea,  with  sufficient  accuracy  and  computational  efficiency. 
Other  key  attributes  desired  in  the  numerical  model  included,  but  were 
not  limited  to: 

a)  flexibility  in  handling  arbitrary  structure  shapes  and  water 
depths  ; 

b)  availability  of  the  code  for  unrestricted  use  by  the  Corps  of 
Engineers;  and 

c)  a  track  record  of  acceptable  performance. 

In  the  case  that  a  model  meeting  all  these  requirements  could  not  be 
found,  the  study  would  recommend  models  which  could  be  modified  and/or 
improved  upon  to  meet  the  Corps  needs. 


2 .  Existing  Models 

An  initial  review  of  the  current  literature  disclosed  at  least  24 
separate  models  dealing  with  the  hydrodynamic  processes  involving  a 
floating  body  in  the  ocean  environment.  These  are  listed  in  Table  1. 
Although  further  searches  will  undoubtedly  expand  the  list,  it  became 
apparent  that  the  24  models  could  be  grouped  into  the  five  principal 
categories  listed  below,  depending  upon  the  approaches  used  in  evaluating 
the  velocity  potential: 
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TABU?  1 

LIST  OF  NUMERICAL  MODELS  FOR  ANALYSIS  OF  FLOATING  BREAKWATER 


NUMERICAL 

APPROACHES  AUTHORS 

FLOATING  STWJCTORE 

APPLICABLE 

WATER 

DEPTH 

SHAPE 

CR0SS- 

SBCTIGN 

DD®*- 

SION 

Eigen 

1  Mei  and  Black  (1969) 

— 

□ 

2 

Constant 

Function 

2  Garrett  ( 1971 ) 

1 

0 

3 

M 

Expansion 

3  Black  et  al.  (1971) 

-  1 

o 

2,3 

II 

Method 

*4  Ijiim  et  al.  (1972a) 

— 

□ 

2 

II 

(EFhM) 

5  Ijima  et  al.  (1972b) 

1 

0 

3 

II 

*6  I to  and  Chiba  (1972) 

— 

a 

2 

II 

7  Yeung  ( 1981 ) 

1 

0 

3 

II 

Green 

8  Kim  (1965) 

A 

3 

Constant 

Function 

*9  Garrison  (1974,  ’78) 

A 

3 

II 

Method 

*10  Adee  ard  Martin  (1974) 

— 

A 

2 

Infinite 

(GFM) 

*11  Georgiadis  and  Hartz 

— 

A 

3 

Variable 

(1982) 

Hybrid 

*12  Ijima  et  al.  (1975) 

1 

a  ° 

3 

Constant 

Green 

*13  Ijima  et  al.  (1976) 

— 

A 

2 

Variable 

Function 

14  Yamamoto  et  al.  (1980) 

— 

A 

2 

II 

Method 

*15  Kiyokava  et  al.  (1983) 

— 

A 

2 

•1 

(HGFM) 

Finite 

16  Zienkiewicz  et  al. 

A 

3 

Variable 

Element 

(1969,  ’77) 

Method 

17  Chen  et  al.  (1974) 

1 

A 

2 

ll 

(FEM) 

*18  Yue  et  al.  (1978) 

A 

3 

li 

19  Hara"Tl978) 

A 

3 

II 

*20  Bando  et  al.  (1984  a,  b) 

A 

3 

ll 

OTHERS 

21  Ursell  (1949) 

_ 

0 

2 

Infinite 

22  Yu  and  Ursell  (1961) 

— 

0 

2 

Constant 

23  Wang  (1966) 

Sphere 

3 

ll 

24  Vugts  (1966) 

Various 

2 

Infinite 

Note:  —  =  Horizontal  O  =  Circular 

I  =  Vertical  □  =  Rectangular 

A  =  Arbitrary  <=>  =  Elliptical 


*  Models  retained  for  additional  detailed  evaluation  in  this 
study. 


v.vv.v 
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o  Eigen  Function  Expansion  Method  (EFEM) 
o  Green  Function  Method  (GFW) 

o  Hybrid  Green  Function  Method  (HGFW) 

o  Finite  Element  Method  (FEM) 

o  Others 

The  Eigen  Function  Expansion  Method  (EFEM)  determines  the  velocity 
potential  through  analytical  expressions  in  the  form  of  Eigen  Function 
expansion.  Applicable  only  to  structure  shapes  of  simple  geometry  and 
constant  water  depths,  the  EFEM  is  the  simplest  numerical  model  of  the 
five  categories  mentioned.  It  is  also  basically  a  two-dimensional  model. 

The  Green  Function  Method  (GFM)  arrives  at  the  velocity  potential  by 
solving  an  integral  eguation  which  contains  the  Green  Function.  The 
integral  equation  is  defined  on  the  structure  surface  from  which  the 
terms  "boundary  integral"  or  "boundary  element"  methods  derive.  Since 
the  numerical  evaluation  of  the  governing  integral  equation  is  performed 
on  the  discretized  elements  of  the  structure  surface,  the  computational 
load  in  the  Green  Function  Method  could  become  quite  burdensome  depending 
upon  the  complexity  of  the  structure  shape. 

The  Hybrid  Green  Function  Method  (HGFM)  combines  some  features  of 
both  the  Green  Function  Method  and  the  Eigen  Function  Expansion  Method. 
The  purpose  is  to  reduce  the  computational  burden  in  cases  of  complex 
breakwater  geometry  and  variable  water  depths.  In  this  method,  the  domain 
of  interest  is  separated  into  interior  and  exterior  sub-regions  to  sim¬ 
plify  the  boundary  conditions,  and  the  velocity  potential  is  expressed  in 
each  sub-region  using  both  the  Eigen  Function  series  and  a  Green 
Function. 


The  Finite  Element  Method  (FEM)  determines  the  velocity  potential  by 
evaluating  the  governing  integral  equation  in  the  entire  domain  of  the 
fluid.  In  the  Finite  Element  Method,  therefore,  the  entire  domain  is 
discretized  into  small  elements,  while  in  the  Green  Function  Method  only 
the  structure  boundary  is  discretized.  Many  existing  FEM  models  employ 
some  degree  of  hybrid  technique  in  order  to  relax  the  radiation  condition 
imposed  at  infinity  and  reduce  the  computational  load.  This  modifica¬ 
tion,  sometimes  called  the  "Hybrid  Element"  method,  could  confine  the 
domain  to  be  discretized  in  a  limited  region  by  employing  analytical 
representations  in  the  far  field. 

Other  Models  typically  employ  a  series  of  singular  functions  to 
represent  the  velocity  potentials  associated  with  body  motions.  This 
method  was  the  forerunner  among  various  approaches  which  were  used  to 
obtain  the  added  mass  and  damping  coefficients  in  ship  motion  problems. 

3.  Mathenatical  Procedures 

3 . 1  TVo  Assumptions 

In  describing  wave  motion  in  the  presence  of  a  floating  break¬ 
water,  there  are  certain  mathematical  features  common  to  all  the 
numerical  models.  One  is  the  use  of  velocity  potential?  another  is 
the  assumption  of  linearity  for  fluid  motion. 

Ti  e  use  of  velocity  potential  is  predicated  on  the  assumption 
that  the  fluid  motion  is  irrotational.  The  assumption  guarantees 
the  existence  of  a  scalar  function  called  "velocity  potential": 

1  (x,  y,  z)  ( 1 ) 


which  satisfies  *h'-  following  condition: 


Thus,  once  the  velocity  potential  <P  is  determined,  all  three 
unknowns  u,  v  and  w  can  be  determined  readily  through  Equation  (2). 
This  is  convenient:  instead  of  pursuing  three  unknowns  simultane¬ 
ously,  one  may  seek  to  determine  a  single  unknown,  <t>  .  In  order  to 
determine  the  velocity  potential,  one  must  solve  a  continuity  equa¬ 
tion,  which,  in  the  irrotational  fluid  motion,  takes  the  form  of 
Laplacian: 


V 


=  o 


(3) 


under  proper  boundary  conditions.  It  is  mainly  in  the  procedures 
for  solving  Eq.  (3)  that  the  major  distinctions  between  different 
numerical  models  arise.  The  remaining  unknown,  pressure  (p),  can 
now  be  determined  by  substituting  <p  into  the  general  Bernoulli 
f  ^ —at ion  (i.e.,  pressure  equation)  of  the  form: 


grad 


+ 


0  (4) 


in  which 


q2  -  u2  +  v2  +  w2  =  o  (linearity  assunption) 
n  =  potential  of  the  system 


p  =  pressure 
P  =  fluid  density 


_  3  ,3  ,3 

3x  3y  3z 


With  velocity  potentials  and  pressure  thus  determined,  the  equations 
of  motion  are  completely  described,  which  can  now  be  solved  for  body 
response  and  then  for  both  wave  height  distribution  and  forces  on 
the  structure.  Although  there  are  some  variations,  most  of  the 
existing  numerical  models  follow  the  computational  steps  outlined 
above. 

Another  feature  which  is  common  to  the  various  numerical  models 
is  the  assumption  that  the  fluid  system  is  linear.  Physically,  this 
assumption  means  that  the  amplitudes  of  the  wave  and  structural 
motions  are  small.  Mathematically,  this  assumption  means  that  the 
velocity  potentials  are  additive.  For  instance,  if  x  denotes 
the  velocity  potential  associated  with  a  unit  heave  of  the  struc¬ 
ture,  the  velocity  potential  associated  with  an  arbitrary  heave  Xj 
will  be  xt  <p1  .If  we  let  $2  denote  the  velocity  potential 
associated  with  a  unit  sway  of  the  structure,  the  velocity  potential 
arising  from  simultaneous  heave  Xj  and  sway  x2  will  be 
X]  £  1  +  x2  $2  «  and  so  on. 

f'sing  the  additive  property  of  the  velocity  potential  in  a 
linear  system,  one  may  first  derive  velocity  potentials  for  all  unit 
motions  and  then  combine  them,  linearly,  as  needed  to  match  the  real 
scenarios.  There  are  usually  eight  such  velocity  potentials: 
f0  for  incident  wave,  throuqh  ,  respectively,  for 


the  heave,  sway,  surge,  roll,  pitch  and  yaw  motions  of  the  structure 
in  calm  water,  and  <p7  for  scattered  wave.  While  the  linearity 
assumption  basically  implies  that  the  motions  involved  are  small, 
the  solutions  are  usually  able  to  describe  moderate-sized  motions 
adequately  as  well. 

3.2  Eigen  Function  Expansion  Method  (EFEM) 


The  EFEM  can  be  considered  the  simplest  numerical  model  among 
the  various  models  discussed  in  this  report,  and  thus  is  applied 
only  to  simple  shaped  structures  in  constant  water  depth.  The  flow 
diagram  for  the  EFEM  is  shown  in  Figure  1.  This  method  seeks  to 
determine  the  velocity  potential  in  closed  form  by  considering 
simplified  boundary  conditions.  For  this  purpose,  the  domain  of 
interest  is  divided  into  two  sub-regions  of  simple  geometry  as  shown 
in  Figure  2:  the  "interior"  sub- region  beneath  the  floating  body 
bounded  by  the  body,  the  sea  bottom  and  the  artificial  borders,  and 
the  "exterior"  sub-region  with  the  free  surface  extending  to  infin- 


In  both  of  the  sub-regions,  the  velocity  potential  can  be 
represented  by  the  form  of  eigen  function  expansion  with  analytical 
expressions  which  satisfy  respective  boundary  conditions  except  at 
the  artificial  boundary.  These  expansions  include  an  infinite 
number  of  unknown  coefficients  which  are  determined  by  requiring 
continuity  of  pressure  and  velocity  across  the  artificial  boundary 
using  either  the  variational  method  (Mei  and  Black,  1969;  Black  et 
al.,  1971)  or  the  Galerkin  method  (Garrett,  1971;  Ijima  et  al., 
1972a,  b). 


The  EFEM  was  applied  to  a  fixed  floating  body  by  Mei  and  Black 
(1969)  and  Garrett  (1971)  and  to  an  oscillating  floating  body  in  a 
calm  sea  by  Black  et  al.,  (1971)  and  Yeung  (1981).  For  the  case  of 
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j  INTERIOR  SUBREGION  EXTERIOR  SUBREGION 


Figure  1.  Flow  Diagram  for  Eigen  Function  Expansion  Method  (EFEM) 


a  floating  body  in  the  presence  of  waves,  Ijima  et  aU,  (1972a,  b) 
developed  a  model  in  which  velocity  potential  and  body  response  are 
determined  simultaneously. 

The  model  developed  by  Ito  and  Chiba  (1972)  essentially  repre¬ 
sents  an  approximate  method  employing  only  the  first  term  of  an 
infinite  series  in  the  Ijima  model.  Their  method  demonstrated 
fairly  consistent  agreement  with  empirical  experiments. 

3.3  Green  Function  Method  (GFW) 


The  Green  Function,  expressed  G(x,  £  ),  defines  a  velocity 
potential  at  a  position  x  when  a  unit  charge  is  placed  at  a  source 
point  £  .  If  source  strength  is  f(  £  )  instead  of  unity,  velocity 
potential  at  x  should  be  f(  £  )  •  G(x,  £  ).  If  charges  of  strength 
f(  £  )  are  distributed  over  a  boundary  S,  velocity  potential  41  (x)  at 
x  is  expressed  by 


=  J  f(  £  )  •  G(x,  £  )ds 

*/S 


Various  forms  of  Green  function  are  known.  One  may  be  selected 
that  allows  velocity  potential  in  accordance  with  the  above  defini¬ 
tion  to  satisfy  the  continuity  equation,  independent  of  f(  £  ). 


^2  $  m  0 


The  selected  Green  Function  may  or  may  not  satisfy  the  following 
series  of  boundary  conditions,  such  as  boundary  condition  at  sea 
surface: 


S  $ 
dz 


(  w2  /g)  ♦ 


(7) 


boundary  condition  at  sea  bed: 


3  $ 
dn 


0  at  z  -  -h 


(8) 


and  the  radiation  condition  at  r  -*■  »  .  The  velocity  potential 
thus  defined  contains  ambiguity  because  f  (  £  )  is  not  known.  This 
function  is  determined  numerically  by  satisfying  the  remaining 
boundary  conditions  and  by  substituting  the  kinematic  boundary  con¬ 
dition. 


3$ 

3n 


(9) 


into  the  integral  equation  (5).  These  two-step  procedures  complete¬ 
ly  define  the  velocity  potential  $  . 

Kim  (1965)  and  Garrison  (1974,  1978)  presented  fully  three- 
dimensional  Green  function  models  applicable  to  floating  breakwaters 
of  arbitrary  geometry  in  constant  water  depth.  Georgiadis  (1982) 
studied  the  response  of  an  infinite  floating  bridge  subjected  to 
oblique  incident  waves  as  a  two-dimensional  case  using  the  cosine 
projection  for  obliquely  incident  waves.  This  approach  is  unable  to 
handle  wave  diffraction  around  the  edge  of  the  breakwater. 


The  computational  procedures  used  in  the  GFM  are  illustrated  in 
the  flow  diagrams  shown  in  Figures  3,  4,  and  5.  In  Figure  3,  the 
entire  computation  is  categorized  into  two  successive  sequences 
consisting  of  "unit  response  analysis"  and  "body  response  analysis". 

The  unit  response  analysis  determines  the  velocity  potentials 
associated  with  unit  motions  of  the  floating  body  in  calm  water  and 
unit  incident  wave.  The  end  products  of  this  sequence  of  computa¬ 
tions  are  velocity  potentials  for  the  unit  of  heave,  sway,  surge, 
roll,  pitch,  yaw  (  4>  x  through  4>s  ),  velocity  potential  for  scattered 
wave  (  $7  ),  plus  added  mass,  damping  coefficient,  and  exciting 
force.  Detailed  steps  in  the  unit  response  analysis  are  shown  in 
Figure  4. 

The  body  response  analysis  essentially  involves  solving  equa¬ 
tions  of  motion  using  the  output  from  the  preceding  unit  response 
analysis.  As  shown  in  detail  in  Figure  5,  it  first  establishes  the 
finite  motions  of  body  response  x  j  through  x6  (corresponding  to 
heave,  sway,  surge,  roll,  pitch  and  yaw)  under  given  finite  wave 
amplitude  H.  The  final  velocity  potential  is  then  computed  through 
an  additive  procedure  allowed  by  the  linearity  assumption,  such  that 


*  “  E  xi  +  H(*o  +  *7) 


The  remaining  steps  leading  to  wave  height  distribution  and 
wave  forces  on  the  structure  are  straightforward,  as  shown  in  Figure 


A  detailed  description  of  the  typical  Green  Function  Method, 
Garrison  Model,  is  presented  in  Appendix  A. 


UNIT  RESPONSE  ANALYSIS 


PRODUCTS: 

■  VELOCITY  POTENTIALS 


FOR  UNIT  BODY  MOTIONS 

01  ....  06 

FOR  UNIT  INCIDENT  WAVE 

07 

■  ADDED  MASS 

■  DAMPING  COEFF. 

■  EXCITING  FORCE 

~T 

BODY  RESPONSE  ANALYSIS 

_  I  z 

PRODUCTS: 

■  WAVE  HEIGHT  DISTRIBUTION 

■  WAVE  FORCES  ON  STRUCTURE 

■  MOORING  FORCE 


Figure  3.  Flow  of  Computational  Sequences  in  the  Green  Function 
Method  (GFM)  and  Finite  Element  Method  (FF.M) 
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UNIT  RESPONSE  ANALYSIS 


0  3  SURGE 
0„  ROLL 
(J)^  PITCH 

C$)6  Y  AW 

0  ,  SCATTERED  WAVE 


Figure  U.  Flow  Diagram  for  Unit  Response  Analysis  in  Green  Function 
Method  (GFM)  and  Finite  F.lement  Method  (FEM) 
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Figure  5.  Flow  Diagram  for  Body  Response  Analysis 


3.4 


Hybrid  Green  Function  Method  (HGFM) 

This  method  combines  partial  features  of  EFEM  and  GFM,  which 
facilitates  computations  dealing  with  complex  breakwater  geometry 
and  variable  water  depths.  In  the  HGFM,  the  domain  of  interest  is 
partitioned  into  interior  and  exterior  sub-regions,  as  in  EFEM.  The 
velocity  potential  in  each  sub-region  is  then  expressed  using  both 
an  Eigen  Function  series  and  a  Green  Function.  The  flow  diagram  in 
the  HGFM  is  shown  in  Figure  6.  In  comparison  with  the  flow  diagram 
for  the  EFEM  (Figure  1),  it  is  clear  that  the  flows  of  computational 
procedures  in  both  methods  cure  essentially  the  same. 

Credit  for  the  initial  development  of  the  hybrid  approach 
belongs  to  I jima  and  his  associates.  Their  models  handle  three- 
dimensional  cases  with  constant  water  depth  (I jima  et  al.,  1975)  and 
two-dimensional  cases  with  variable  water  depths  (I jima  et  al., 
1976).  The  hybrid  model  introduced  by  Kiyokawa  et  al_.,  (1983) 
essentially  boils  down  to  a  minor  modification  of  the  I  jima  two- 
dimensional  model  allowing  for  easy  determination  of  hydrodynamic 
coefficients  (i.e.,  added  mass,  damping  coefficient,  and  exciting 
force).  A  more  detailed  description  of  the  Ijima  model  (1975)  is 
presented  in  Appendix  A. 

3 . 5  Finite  Element  Method  (FEM) 


The  flow  of  computational  steps  in  the  FEM  is  exactly  the  same 
as  for  the  GFM,  as  shown  in  Figures  3,  4,  and  5.  The  principal 
difference  between  the  two  models  arises  from  the  methods  for  deter¬ 
mining  velocity  potentials  under  unit  excitation  in  the  unit  re¬ 
sponse  analysis.  In  the  FEM,  the  continuity  equation  (  7 2  <j>  =  0  )  is 
transformed  into  integral  equations  which  are  defined  over  the 
entire  domain  of  the  fluid,  in  contrast  to  the  GFM  in  which  the 
governing  integral  equations  are  defined  only  on  the  boundary.  In 
the  FFM,  therefore,  the  entire  domain  is  discretized  into  small 
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Figure  b.  Flow  Diagram  lor  Hybrid  Hreen  Function  Method  iH(-t-M) 


19 


elements,  while  the  discretization  in  the  GFM  applies  only  to  the 
boundary,  as  shown  in  the  examples  in  Figure  7. 

The  initial  FEM  model  proposed  by  Zienkiewicz  (1969)  was 
designed  to  compute  scattered  waves  around  a  fixed  structure  only. 
Hara  (1978)  expanded  this  model  to  allow  for  the  analysis  of 
floating  bodies  including  both  wave  height  distribution  and  forces 
on  the  structure.  Bando  et  al. ,  (1984a,  b)  introduced  further 
improvements  by  incorporating  computational  accuracy  and  efficiency. 
Chen  et  al. ,  (1974)  introduced  a  simplified  version  of  the  FEM  with  a 
long  wave  assumption,  and  Yue  et  al . ,  (1978)  extended  the  long-wave 
model  to  three-dimensional  problems. 

The  accuracy  and  computing  time  in  the  FEM  are  sensitive  to  the 
way  the  radiation  condition  at  infinity  is  treated.  In  order  to 
formulate  the  radiation  condition  into  the  FEM,  the  behavior  of  the 
velocity  potential  in  the  far  field  should  be  expressed  by  a  certain 
analytical  form,  since  the  entire  fluid  domain  extending  to  infinity 
could  not  be  covered  by  usual  finite  elements.  The  Yue  et  al., 
model  uses  an  Eigen  Function  expansion  for  far  field  behavior,  while 
the  Bando  et  al.,  mouel  adepts  a  simple  decaying  function.  In  this 
sense,  both  approximate  methods  may  be  termed  the  Hybrid  Finite 
Element  Method.  In  fact,  Yue  et  al.,  (1978)  call  their  model  the 
"Hybrid-Element"  method. 

3.6  Other  Methods 

In  all  models  of  this  category  the  potential  is  written  as  the 
sum  of  wave  potentials  with  a  singularity  each  satisfying  all  condi¬ 
tions  except  the  body  surface  boundary  condition.  This 
boundary  condition  is  then  used  to  determine  the  unknown 
coefficients  included  in  each  potential.  Since  the  mathematical 
derivation  and  numerical  manipulation  are  complicated,  it  appears 
difficult  to  generalize  this  approach  for  scatterinq  by  a  lonq 


A.  DISTRIBUTION  OF  ELEMENTS  IN  DISCRETIZED  COMPUTATIONAL  DOMAIN 
FOR  FLOATING  CIRCULAR  CYLINDER  IN  FEM 


B.  DISCRETIZED  BOUNDARY  OF  FLOATING  CYLINDER  IN  GFM 


cylinder  in  obliquely  incident  waves  or  for  arbitrary  three-dimen¬ 
sional  hollies. 

This  method  was  intiallv  developed  by  Ursell  (1949),  for  the 
case  of  a  half-submerged  circular  cylinder  in  calm  water  of  infinite 
depth,  then  extended  by  others  to  various  shapes  and  finite  water 
depths.  Wang  (1966)  applied  the  approach  to  heaving  oscillations  of 
a  half  submerged  sphere  in  a  finite  water  depth. 

Randan  Wave  Input 

4. 1  Discretization 

The  directional  wave  spectrum  S(f,  )  is  expressed  as: 


S(f,  0  )  =  S(f)  •  G(f,  ) 


in  which  S(f)  is  the  frequency  spectrum  and  G(f,  v  )  is  the 
directional  function,  so  that 


sm  »  J_ 


S(f,  )d 


and  thus 


J  G(f,  )< 


1  --  1 


.’•Vr  A 


For  the  purpose  of  numerical  computation,  Eq.  (11)  is  dis¬ 
cretized,  such  that 


S(f.  0  )  = 


«  n 

E  I"  A  S(f  ±)  (E  G<fi«  9  j  )  >1 

i  =  1  L  j  =  1  J 


From  the  numerical  point  of  view,  decisions  must  be  made  on  the 
optimal  number  (n)  of  discrete  bands  in  the  directional  variability 
and  that  of  the  frequency  domain  (N).  For  instance,  as  the  number 
of  wave  direction  bands  (n)  is  increased,  the  accuracy  of  the  compu¬ 
tation  will  improve,  but  the  rate  of  improvement  would  diminish  as 
the  number  n  exceeds  some  critical  value.  Thus,  there  exists  an 
optimum  value  for  n  where  the  accuracy  and  computational  economy  are 
best  balanced.  The  same  argument  should  apply  to  the  optimal  number 
of  frequency  bands,  N. 

4.2  Optimum  Discretization  of  Wave  Direction 


Optimum  discretization  of  wave  direction  has  been  investigated 
in  a  series  of  papers  by  Goda  et  al.,  (1975,  'll),  Nagai  (1972  a,  b, 
c),  Nagai  et  al.,  (1974),  and  Takayama  et  al.,  (1978).  These 
Japanese  studies  typically  used  the  Mitsuyasu  directional  function 
(1981)  with  modifications  introduced  by  Goda  and  Suzuki  (1975). 

The  Mitsuyasu  directional  function  contains  the  so-called  wave 
spreading  parameter  S  which  relates  the  spreading  characteristics  to 
wave  frequency.  For  instance,  in  any  sea  state,  wave  spreading  is 
at  its  minimum  near  the  peak  frequency  fp.  It  is  less  in  swell  than 
in  wind  waves,  and  is  progressively  less  with  decay  distance. 
Figure  8  shows  the  directional  function  by  Mitsuyasu  (1981). 
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The  Mitsuyasu  directional  function  is  written  as: 


G(f,  G  )  = 


I  22s-1  (S! )2 

71  (2s): 


2S  fi 

cos  — 
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(15) 
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in  which 


S  = 


-2.5 


max 


(_r) 


'  s  (-f)5 

max  I  p 


for  f  >  ft 


for  f  is  f. 


(16) 
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according  to  Goda  and  Suzuki  (1975).  The  term  Smax  is  the  spreading 
parameter  for  the  peak  frequency  f  =  fp  ,  and  is  expressed  as 


-5  max  =  11.5  (2v  fpU10/g)  2‘5 


(18) 


in  which  IJ  id  is  the  wind  speed  it  a  10-m  height  above  water. 
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In  practical  cases,  Smax  assumes  the  following  values  (Goda  and 
Suzuki,  1975;  Goda,  1985): 


S  max 

n 

H-* 

o 

Wind  wave 

c 

Jmax 

=  25 

Swell  with  a  medium  steepness  after  a  short 

decay  distance 

c 

°  max 

=  75 

Swell  with  a  lew  steepness  after  a  long  decay 

distance 

Discretization  scheme  must  be  adjusted  to  the  Smax  values 
because  the  spreading  characteristics  vary  greatly  due  to  this 
parameter.  For  instance,  when  Smax  is  large,  cane  should  not  be  con¬ 
cerned  over  the  direction  near  ±  */ 2.  Takayama  et  al.,  (1978) 
showed  limiting  wave  directions  0max  in  the  numerical  computation 
to  be  considered  as  the  function  of  Smax  as  follows: 


Another  well-known  directional  function  is  the  SWOP  (Stereo 


Wave  Observation  Program;  see  Kinsman,  1965)  model.  The  wave 
spreading  behavior  in  the  SWOP  model  is  frequency-dependent  and  is 
the  function  of  a  parameter  related  to  the  wind  speed  at  a  5-m 
elevation,  an  ambiguous  parameter.  This  ambiguity  must  be  resolved 
before  the  SWOP  model  can  be  used  for  practical  application. 


In  addition  to  a  spreading  parameter  ^„ax ,  there  exists  another 
parameter  which  influences  the  discretization  of  wave  direction; 


the  location  of  interest  relative  to  the  breakwater.  This  addition¬ 


al  parameter  is  essentially  a  distance  from  the  tip  of  the  breakwat- 


in  which  R  is  the  distance  and  L  the  local  wave  length. 


Takayama  et  al_.,  (1978)  have  shown  that  for  any  value  of  S^ax , 
the  larger  the  parameter  R/L,  the  larger  must  be  the  number  of 
directional  discretization  in  order  to  preserve  an  acceptable  level 
of  accuracy  in  the  computation.  In  the  case  of  a  semi-infinite 
breakwater,  these  investigators  proposed  optimum  values  for  n  as  a 
function  of  Smax  and  R/L  as  shown  in  Table  2.  Their  results 


indicate  that  the  number  of  directional  discretization  can  be  small 


if  the  location  of  interest  lies  close  to  the  tip  of  the  breakwater. 
However,  for  locations  far  from  the  breakwater  tip,  a  greater  number 
of  discretization  is  necessary. 
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TABLE  2.  OPTIMUM  NUMBER  OP  DIRECTIONAL  DISCRETIZATION 
FOR  WAVES  INCTDEKT  TO  SEMI -INFINITE  BREAKWOTR 


R/L 

Snnx 

10 

25 

75 

1000 

5 

8 

8 

6 

6 

10 

18 

15 

12 

9 

20 

20 

19 

16 

10 

50 

25 

23 

20 

12 

80 

30 

25 

23 

13 

Takayama  et  al.(  (1978)  suggested  a  formula  giving  the  optimum 
value  of  n  as  the  function  of  Smax  and  R/L  as  follows: 


n  =  {  19.1  -  3.8  log  (Smax  )  {  log  (R/L) 


Figure  9  shows  an  example  of  computation  of  diffraction  coeffi¬ 
cients  behind  a  semi-infinite  breakwater  for  three  directional 
bands:  5°  (n  =  36),  10°  (n  =  18),  and  20°  (n  =  9).  The  principal 

direction  of  the  wave  was  perpendicular  to  the  breakwater.  In  this 
case,  the  directional  discretization  was  obviously  inadequate  at  n  = 
9  but  too  fine  at  n  =  36.  Using  Eq.  (19),  the  optimum  value  for  n 
is  determined  as  n  =  29.11  ~  30  (i.e.,  6°)  for,  say,  R/L  =  80.  For 
R/L  =  40  (i.e.  a  position  closer  to  the  tip  of  the  breakwater),  one 
gets  n  =  24.5  ~  25  (i.e.  7C). 


B  MONOCHROMATIC  WAVE 


A  RANDOM  WAVE  WITH  MITSUYASU 
DIRECTIONAL  FUNCTION 

Takayama  et  al.  (1978) 


Figure  9.  Contrasting  Diffraction  Patterns  Due  to  Ditterent  Input  Waves 
on  Semi-Infinite  Breakwater:  A.  Random  Wave  With  Directional 
Spreading.  B.  Monochromatic  Wave. 


4.3  Discret izatian  of  Fr< 


equency  Bands 

Takayama  et  al.  (1978)  performed  a  numerical  experiment  in 
which  fluctuations  in  the  diffraction  coefficient  Kj  were  investi¬ 
gated  as  a  function  of  the  number  of  discretized  frequency  bands,  N, 
in  the  case  of  a  semi-infinite  breakwater.  The  experiment,  in  which 
the  directional  discretization  was  held  at  n  =  36,  showed  that  the 
diffraction  coefficients  K<j  remained  essentially  stable  for  values 
of  N  greater  than  2.  The  diffraction  coefficient  was  also  found  to 
fluctuate  little  with  variations  in  S  max  and  R/L.  This  meant  that 
the  computational  accuracy  was  far  less  dependent  on  the  discre¬ 
tization  of  frequency  bands  than  that  of  the  wave  direction.  In 
this  experiment,  for  instance,  if  the  computational  stability 
was  defined  as 


Kd  -  K<T 


G 


(20) 


in  which  K<j  and  K<j'  are  diffraction  coefficients  for  discretization 
N  and  N'  (N  <  N'),  respectively,  N  was  as  low  as  N  =  5  to  ensure 

e  <  2%  (21) 


Furthermore,  Nagai  (1972)  recommends  that  the  discretization  of 
frequency  bands  must  satisfy  the  following  two  conditions  (Figure 
10). 


1°  The  non-dimensional  spectral  energy  contained  in  each 
discretized  band  should  be  the  same: 


1  r 

A  E  -  —  J  S(f*)  df *  (22) 

0 

30 


in  which  f*n  is  the  non-dimensional  cutoff  frequency  and  f  =  f/fp. 

2°  The  representative  frequency  fc*  in  each  discretized  fre¬ 
quency  band  should  be  determined  in  such  a  way  that  the 
second  moment  in  each  band  around  fc*  is  the  same  as  the 
second  moment  relative  to  f  =  0,  i.  e., 


f 


*2 

c 


AE 


f*2S(f*)df* 


(23) 


in  which 


Nagai  (1972)  derived  analytical  expressions  for  AE  and  fc*  in  the 
case  of  Pierscn-Mosoowitz  and  Bretschneider  spectra  ( see  Figure  10 ) . 
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In  addition  to  these  criteria,  which  lead  to  smaller  dis¬ 
cretized  bands  around  the  peak  frequency,  attention  should  also  be 
paid  to  the  resonance  frequency  of  a  floating  body.  The  body  re¬ 
sponse  could  grow  sharply  at  the  resonance  frequency  even  though  the 
incident  wave  energy  at  this  frequency  may  be  relatively  small. 
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Figure  10.  ExampLe  of  Discretization  of  Frequency  Bands  (Nagai,  1972) 
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Evaluation  of  Models 


5.1  Initia 


The  comprehensive  list  shown  in  Table  1  was  prepared  to  include 
numerical  models  dealing  with  floating  bodies  in  general.  Some  of 
these  models  are  not  exactly  applicable  to  floating  breakwaters 
without  substantial  modifications  and  hence  should  be  removed  from 
further  consideration  in  this  study. 

In  the  group  of  Eigen  Function  Expansion  Method,  all  but  Ijima 
et  al^.,  (1972a)  and  Ito  and  Chiba  (1972)  were  removed.  Models  by 
Mei  and  Black  (1969)  and  Garrett  (1971)  apply  to  a  fixed  floating 
body  and  those  by  Black  et  al.,  (1971)  and  Yeung  (1981),  only  to  an 
oscillating  floating  body  in  a  calm  sea.  Another  model,  Ijima  et 
al.,  (1972b)  treats  a  vertical  cylinder  which  obviously  departs  from 
the  usual  geometry  of  the  floating  breakwater.  The  model  by  Ijima 
et  al.,  (1972a),  which  is  retained  in  this  group,  deals  with  an 
elongated  floating  structure  and  hence  represents  a  two-dimensional 
case  of  a  floating  breakwater  with  constant  water  depth.  The  model 
by  Ito  and  Chiba  (1972)  is  an  abridged  version  of  Ijima  et  al., 
(1972a),  in  which  only  the  first  terms  of  the  infinite  Eigen  Func¬ 
tion  series  were  retained  to  achieve  computational  economy. 


In  the  group  of  Green  Function  Method,  there  exist  a  number  of 
three-dimensional  Green  Function  models  similar  to  those  introduced 
by  Kim  (1965)  and  Garrison  (1974,  *78).  But  we  decided  to  retain 
only  the  Garrison  model  inasmuch  as  the  differences  among  them  are 
small.  The  model  by  Georgiadis  and  Hartz  (1982)  was  retained  to 
enhance  evaluation,  because  this  model  has  the  added  capability  of 
incorporating  short -crested  waves  as  input. 
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In  the  group  of  Hybrid  Green  Function  Model,  the  model  de¬ 
scribed  by  Yamamoto  et  al.,  (1980)  is  exactly  the  same  as  that  of 
Ijima  et  aL,  (1976),  and  was  therefore  removed. 

In  the  group  of  Finite  Element  Model,  the  models  by  Zienkiewicz 
et  al.,  (1969,  '77)  and  Hara  (1978)  were  removed,  because  they  are 
superceded  by  the  more  updated  version  described  by  Bando  et  al., 
(1984a,  b).  All  the  models  in  the  category  of  "Others"  were  removed 
from  further  evaluation,  because  these  models  were  intended  for 
analysis  of  ship  motions  in  the  seaway,  and  were  not  completely 
applicable  to  a  floating  breakwater. 

As  a  result  of  the  screening  process  described  above,  the 
following  ten  models  were  retained  for  further  evaluation: 

Eigen  Function  Expansion  Models: 

Ijima  et  al_. ,  (1972a) 

I to  and  Chiba  (1972) 

Green  Function  Models: 

Garrison  (1974,  '78) 

Adee  and  Martin  (1974) 

Georgiadis  and  Hartz  (1982) 

Hybrid  Green  Function  Models: 

Ijima  et  al.,  (1975) 

Ijima  et  al.,  (1976) 

Kiyokawa  et  al.,  (1983) 


Finite  Element  Models: 


Yueetal. ,  (1978) 

Bando  et  al . ,  ( 1984a, b) 

5.2  Criteria  for  Evaluation 

The  models  were  evaluated  against  nine  criteria: 


.  Accuracy  (15%) 

C.P.U.  Time  (15%) 

.  Storage  Requirement  (10%) 

.  Input  Preparation  (15%) 

.  Diffracticn  Effect  (15%) 

.  Directional  Wave  (10%) 

.  Breakwater  Shape  (12%) 

.  Water  Depths  (8%) 

.  Code  Availability  (0%) 


The  numbers  within  the  parentheses  represent  the  percent  weight 
assigned  to  each  criteria. 

We  considered  that  the  computational  performance  of  the  model 
should  be  given  more  than  50%  of  the  overall  significance.  Accord¬ 
ingly,  we  assigned  a  weight  of  55%  to  the  first  four  criteria: 
accuracy,  C.P.U.  (central  processing  unit)  time,  computer  storage 
requirements,  and  the  work  required  to  generate  the  input  data 
(particularly  the  work  required  for  discretization).  The  functional 
attributes  of  the  model,  consisting  of  the  capabilities  to  incorpor¬ 
ate  three-dimensional  effects  (i.e.,  diffracticn  processess  and  di¬ 
rectional  wave  spectrum  for  the  input  wave)  and  the  flexibility  to 
accomodate  various  breakwater  shapes  and  water  depths,  were  assigned 
a  combined  total  weight  of  45%.  The  code  availability  was  not 


assigned  a  weight  because  it  was  safe  to  assume  that  the  code  would 
be  available  in  any  case. 

The  tally  of  the  scores  for  each  screened  model  is  summarized 
in  Table  3. 

5 . 3  Quantitative  Evaluation 

Regarding  the  accuracy  of  computation,  all  the  models  with  the 
exception  of  Ito  and  Chiba  (1972)  and  Georgiadis  and  Hartz  (1982), 
were  given  the  maximum  score  of  15%.  The  model  by  Ito  and  Chiba 
(1972)  employed  only  the  first  term  of  an  infinite  Eigen  Function 
series  for  velocity  potential  in  the  Ijima  et  al.,  (1972a)  model. 
Although  this  was  an  approximate  solution,  they  reported  fairly 
close  agreement  between  numerical  results  and  the  experimental  data. 

In  the  model  by  Georgiadis  and  Hartz  (1984),  the  radiation 
condition  which  is  only  valid  at  infinity  is  prescribed  on  the 
artificial  boundary  at  a  finite  distance  from  the  breakwater.  Al¬ 
though  this  artifact  undoubtedly  provides  computational  expediency, 
the  influence  of  this  approach  on  the  accuracy  of  the  results  is 
unknown  and  should  be  examined. 

The  criteria  of  computing  time  (C.P.U.  time),  computer  storage 
requirements  and  input  data  generation  relate  to  computational  eco¬ 
nomy.  The  three-dimensional  FEM  (Yue  et  al.,  1978;  Bando  et  al., 
1984a,  b)  models  and  the  three-dimensional  Green  Function  models 
(Garrison  1974,  '78)  scored  poorly  with  respect  to  computational 
economy.  This  is  mainly  because  these  models  require  discretization 
in  the  three-dimensional  domain  of  interest,  usually  a  labor-inten¬ 
sive  effort.  Out  of  the  maximum  possible  score  of  40%  for  computa¬ 
tional  economy,  the  FEM  models  by  Yue  et  al.,  (1978)  and  Bando  et 
al.,  (1984a,  b)  were  given  only  12%  each,  and  the  GFM  model  by 
Garrison  (1974,  '78)  only  18%.  The  three-dimensional  Hybrid  Green 


TABLE  B .  QUANTITATIVE  EVALUATION  OF  SCREENED  MODELS 


Function  Model  (Ijima  et_  al.,  1975)  scored  29%,  demons! rat  inq  a 
considerable  advantage  over  the  three-dimensional  FEM  and  GF1M  mod¬ 
els.  In  fact,  two-dimensional  models  scored  30%  or  above  on  compu¬ 
tational  economy  for  obvious  reasons. 

In  Table  3,  all  the  three-dimensional  models,  marked  *,  possess 
the  ability  to  accomodate  directional  wave  spectra  and  diffraction 
effect  around  the  edge  of  the  breakwater.  The  model  by  Georgiadis 
and  Hartz  (1982)  is  not  exactly  a  three-dimensional  model,  since  it 
cannot  handle  diffraction  effect,  although  it  can  accommodate  the 
short-crested  oblique  waves.  This  model  deals  with  an  infinitely 
long  floating  structure,  essentially  a  two-dimensional  approach. 

The  models  using  the  Eigen  Function  Expansion  Method  have 
limitations  in  terms  of  accommodating  versatile  breakwater  geometry. 
The  hybrid  model  by  Ijima  et  al.,  (1975)  dealt  with  cylinders  with 
elliptical  and  rectangular  sections,  which  can  be  elongated  into  a 
horizontal  dimension  to  simulate  the  usual  oblong  shape  of  the 
floating  breakwater.  The  FEM  and  GFM  models  typically  have  the 
ability  to  accomodate  arbitrary  three-dimensional  breakwater  geomet¬ 
ry. 

Whereas  the  FEM  models  can  handle  variable  water  depths,  the 
GFM  models  cannot  because  the  Green  Function  is  usually  derived  for 
a  constant  water  depth.  However,  this  typical  limitation  of  the  GFM 
model  can  be  overcome  through  modifications  such  as  those  incorp¬ 
orated  by  Georgiadis  and  Hartz  (1982),  Ijima  et  al_.,  (1976)  and 
Kiyokawa  et  al.,  (1983).  The  model  by  Adee  and  Martin  (1974)  only 
applies  to  infinite  water  depth,  hence  is  not  appropriate  to  a 
finite  depith  where  tlie  floating  breakwater  is  installed. 


In  order  to  gain  a  further  understanding  in  computational  time, 
a  comparison  of  the  CFU  time  among  the  models  has  been  made.  Infor¬ 
mation  on  the  performance  of  six  programs  based  on  the  following 
numerical  models  was  obtained.  The  results  of  the  Yue  et  al.,  model 
were  found  in  their  paper  (1978).  The  others  were  obtained  through 
private  communication  with  Mr.  Takahashi  of  the  Kajima  Corporation. 

Eigen  Function  Expansion  Models: 

Ijimaetad.,  (1972a) 

I to  and  Chiba  (1972) 

Green  Function  Models: 

Garrison  (1974,  *78) 

Finite  Element  Models: 

Yue  et  al.,  (1978) 

Bandoetal. ,  (1984a,  b) 

To  arrive  at  a  standard  of  comparison,  the  response  of  a  float¬ 
ing  dock  in  a  simple  harmonic  wave  field  was  computed  by  each  pro¬ 
gram  (see  Figure  11).  Table  4  summarizes  the  CPU  time  per  single 
wave  frequency  as  a  function  of  the  number  of  unknowns  solved  by 
each  model.  The  number  of  unknowns  is  essentially  associated  with 
the  number  of  simultaneous  equations  to  be  solved,  and  in  the  case 
of  the  GFM  and  the  FEM,  it  is  further  related  to  the  number  of  dis¬ 
cretized  elements. 


The  computer  used  by  Takahashi  was  a  Hitachi  M-280H,  considered 
equivalent  to  an  IBM  3081K.  The  computer  used  by  Yue  et_  ab,  (1978) 
is  apparently  a  main-frame  variety,  but  its  exact  model  is  unknown. 

Both  programs  based  on  the  EFEM  demonstrated  the  remarkable 
computational  efficiency,  showing  a  faster  CPU  time  by  1  to  3  orders 
of  magnitude.  The  discrepancy  in  CPU  time  between  the  two  GFM 
programs  is  attributed  to  the  different  numerical  evaluation  methods 
for  the  Green  Function.  The  second  GFM  program  is  expected  to 
contain  an  advanced  subroutine  which  speeds  up  the  evaluation  of 
Green  Function. 

Using  the  FEM  model  by  Bando  et  al_.,  computations  for  the 
different  numbers  of  unknowns  resulted  in  widely  varying  CRJ  times 
as  shown.  As  a  rule  of  thumb,  the  CPU  time  in  the  FEM  is  propor¬ 
tional  to  the  square  of  the  number  of  unknowns.  In  this  case,  where 
the  number  of  unknowns  is  reduced  to  a  practical  minimum  of  400,  the 
CPU  time  is  estimated  to  achieve  the  level  comparable  to  that  of  the 


Figure  12  shows  the  body  response  computed  by  the  I  to  and  Chiba 
and  the  Bando  et  al.,  models.  Agreement  between  the  two-dimensional 
and  the  three-dimensional  models  is  fairly  close  in  terms  of  body 
response . 


METHOD 

MODEL  # 

CPU  (sec) 

NO.  OF  UNKNOWNS 

IN  TABLE  3 

EFEM  (Ito)*a 

6 

0.005 

5 

EFEM  (Ijima) 

4 

1 

3  +  ? 

GFM  ( Garrison ) 

9 

100 

45 

GEM*C 

9 

11 

45 

FEM  (Bando) 

20 

1150 

3893 

FEM  (Bando) 

20 

300 

2165 

FEM  (Bando) 

20 

10  -  20*^ 

400 

FEM  (Yue)*e 

18 

200*f 

435 

*d  :  Estimation 

*e  :  Exciting  force  only 

*f  :  Computed  in  1978 


Written  by  Bando 
Written  by  Bando 

Ccnmercial  package  sold  by  Veritec 


5.5  Costs 


Program  cost  is  one  of  the  most  important  parameters  in  the 
selection  of  the  best  numerical  model  for  CERC.  The  following  Table 
summarizes  the  cost  information  obtained  through  communications  with 
Professor  Ijima  of  Kyushu  University,  Professor  Newman  of  MIT,  the 
MIT  Patent  Office,  Dr.  Garrison ,and  the  Kajima  Corporation  in  Japan. 

TABLE  5.  OOSTS  OF  VARIOUS  PROGRAMS 


PROGRAM 

AUTHOR  OR  OWNER 

COST 

DELIVERY  FORM 

EFEM  (6) 

Kajima  Corp. 

$1,000 

Source  Code 

HGFM  (12) 

Ijiira,  T. 

$4, 000-$ 5, 000 

Source  Code 

GFM  (9) 

Garrison,  C.  J. 

$35,000 

Source  Code 

GFM  (9) 

Kajima  Corp. 

$6,000 

Source  Code 

FINGREEU 

Newman,  J.  N. 

$200 

Object  Module 

FINGREEN 

Newman,  J.  N. 

$10,000 

Source  Code 

_ i 

The  number  within  the  par  an  theses  denotes  the  model  number  from 
Table  3. 

The  Eigen  Function  Expansion  Method  (6)  based  on  the  Ito  and 
Chiba  two-dimensional  model,  was  written  by  Bando  of  Kajima.  Al¬ 
though  the  license  fee  for  the  GFM  (9)  program  owned  by  Garrison  is 
$20,000,  a  higher  fee  has  been  set  for  CERC  because  the  program  is 
expected  to  be  used  by  quite  a  large  number  of  people.  Kajima's 
GFM  (9)  program  was  written  by  Bando  after  the  original  paper  by 
Garrison.  Governmental  agencies  can  obtain  the  object  module  of  the 


program  "FINGREEN"  for  a  handling  fee  of  only  $200,  but  must  pay 
$10,000  for  the  source  code. 

Reocnraendaticns 

6.1  Discussion 

The  Hybrid  Green  Function  Model  by  Ijima  et  al.,  (1975)  scored 
well  in  all  the  criteria  tested,  achieving  the  highest  tally  of  84%. 
It  also  outscored  the  second-ranking  model  by  Garrison  (1974,  *78) 

by  a  clear  margin  of  9%.  Comparing  these  two  high-ranking  models, 
we  find  that  the  Ijima  et  al.,  (1975)  model  excels  in  computational 
economy  (29%  versus  18%),  while  the  Garrison  (1974,  '78)  model  is 
somewhat  superior  in  terms  of  the  ability  to  accommodate  breakwater 
shapes  (12%  versus  10%). 

Recently,  an  efficient  algorithm  was  developed  (Newman,  1985) 
to  expedite  the  computational  speed  in  the  evaluation  of  the  Green 
Function.  Packaged  into  a  subroutine  named  "FINGREEN",  this  method 
replaces  the  numerical  evaluation  of  the  relevant  integrals  with 
multidimensional  approximations  in  economized  polynomials.  A  sub¬ 
stantially  faster  computational  speed  is  said  to  be  achieved  by  this 
method  than  by  conventional  methods  based  on  direct  numerical  inte¬ 
gration.  We  have  not  had  the  opportunity  of  testing  whether  or  not 
the  computational  efficiency  of  the  Garrison  model,  when  coupled 
with  "FINGREEN",  would  exceed  that  of  the  Ijima  et  al.,  (1975) 
model.  A  more  detailed  description  of  "FINGREEN"  is  presented  in 
Appendix  B. 

Three  other  models,  those  by  Georgiadis  and  Hartz  (1982),  Yue 
et  al.,  (1978)  and  Bando  et  al.,  (1984a,  b),  scored  an  identical 
72%,  ranking  slightly  behind  the  Garrison  model  which  scored  75%. 
Of  these  three  third-ranking  models,  the  Georgiadis  and  Hartz  (1982) 
modeL  is  not  strictly  three-dimensional,  since  it  cannot  incorporate 
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the  wave  diffraction  effect.  The  models  by  Yue  et  al.,  (1978)  and 
Bando  et  al.,  (1984a,  b)  are  very  similar  to  each  other. 

The  mooring  force  is  one  of  the  most  important  sources  of 
information  in  the  design  of  floating  breakwaters.  Although  this 
force  is  computed  automatically  in  the  model,  the  numerical  result 
is  highly  dependent  cn  the  spring  constant  given  by  the  users.  The 
spring  constant,  which  is  the  gradient  of  the  displacement-force 
curve,  depends  on  the  magnitude  of  its  displacement,  hence  varies 
continuously  throughout  the  motion  of  the  floating  body.  Figure  13 
shows  an  example  of  this  non-linear  relationship  in  the  case  of  a 
single  catenary-shaped  mooring  cable.  A  practical  idea  is  to  pre¬ 
select  a  spring  constant  by  anticipating  the  range  in  which  the 
displacement  would  occur.  Another  problem  is  the  conjugation  of 
mooring  forces  in  a  multiple  cable  system.  To  date,  little  is  known 
as  to  hew  to  distribute  the  total  computed  mooring  force  to  individ¬ 
ual  cables. 


Presently,  state-of-the-art  models  with  floating  breakwater 
applications  are  continuously  being  updated  and  improved.  The 
preceding  discussion  was  based  on  observations  made  during  the 
inmediate  time  this  review  was  performed. 

6.2  Reccnrnendations 

We  conclude  that  the  Hybrid  Green  Function  Model  by  Ijima  et 
al. ,  (1975)  has  the  best  overall  attributes.  It  not  only  achieved 
the  highest  total  score  with  a  clear  margin  over  the  next  best 
model,  it  also  scored  evenly  on  all  the  criteria  tested.  The  Ijima 
et  al. ,  (1975)  model  has  been  developed  as  a  research  tool,  and  is 
not  in  a  user-friendly  form.  To  make  the  program  fit  for  routine 
use,  refinement  of  program  code,  preparation  of  program  documenta¬ 
tion  and  a  quality  control  review  should  be  carried  out.  Prof. 
Ijima  of  Kyushu  University,  the  co-author  of  the  model ,  expressed 
willingness  to  cooperate  in  these  tasks.  The  Ijima  et  al . ,  (1975) 
model  can  be  operated  on  a  minicomputer . 
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For  the  purpose  of  obtaining  a  quick  diagnostic  assessment 
before  committing  a  sophisticated  three-dimensional  analysis,  the 
two-dimensional  model  by  Ito  and  Chiba  deserves  special  attention. 
This  model  can  be  operated  on  a  microccmputei  and  is  extremely 
inexpensive. 

Ongoing  improvements  of  state-of-the-art  numerical  modeling 
tools  require  users  to  seek  information  concerning  the  most  recent 
advancements .  Future  models  will  be  improved  and  related  litera¬ 
ture  will  be  continuously  updated. 

In  summary,  the  following  recommendations  are  made  at  this 

tine: 

1)  Acquire  the  Ijima  et  al. ,  (1975)  model  and  adapt  it  into  a 
user-friendly  form.  The  program  itself  is  available  at  no 
cost,  but  the  tasks  for  necessary  adaption  would  cost 
$4,000  to  $5,000. 

2)  Cbnduct  test  runs  of  the  Ijima  et  al.  model  against 
selected  scenarios  of  prototype  tests  by  CERC  and  the 
Japanese  Ministry  of  Transport.  In  these  test  runs,  random 
waves  with  directional  variability  shall  be  used  as  input 
waves. 

3)  The  two-dimensional  Ito  and  Chiba  model  is  worth  considera¬ 
tion  as  a  tool  for  speedy  diagnostic  assessment. 

4)  If  additional  financial  resources  are  available,  the 
combination  of  Garrison's  GFM  model  and  the  "FINGREEN" 
subroutine  may  be  worth  consideration. 

5)  Efforts  to  improve  the  computation  of  the  mooring  force 
shall  be  made  as  part  of  the  test  runs  by  paying  special 
attention  to  the  existing  field  data. 
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Appendix  A: 


Detailed  Description  of  Selected  Models 

Ijima  et  al.,  (1975)  Model  (Hybrid  Green  Function  Model) 

In  this  method  the  domain  of  interest  is  divided  into  two  sub- 
regions  of  simple  geometry  as  shown  in  Figure  2:  the  "interior"  sub- 
region  beneath  the  floating  body  bounded  by  the  body,  the  sea  bottom  and 
the  artificial  borders,  and  the  "exterior"  sub-region  with  the  free 
surface  extending  to  infinity. 

In  each  of  the  sub-regions,  the  velocity  potential  can  be 
represented  by  the  horizontal  distribution  functions  (Green  Function)  and 
the  Eigen  Functions  in  the  vertical  direction.  The  velocity  potential  in 
the  exterior  sub-region  satisfies  any  conditions  other  than  the  boundary 
conditions  at  infinity  and  at  the  artificial  border,  while  the  velocity 
potential  in  the  interior  sub-region,  including  some  components  of  body 
motion,  satisfies  all  the  boundary  conditions  but  the  one  at  the  artifi¬ 
cial  border.  They  are  written  in  the  following  form: 


0E  (x,y,z) 


=  JJL[{f0 (x,y)  +  (x,y) } 

2  o 


cosh  k(z  +  h) 
cosh  kh 
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cos  kn(z+h)j 
cos  knh 
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kh  tanh(kh)  =  -  knh  tan(knh)  -o2 h/g 
f0 (x,y)  =  -i  exp  {-ik( xcos  w  +  y  sin  w )} 
q'=  1  -  q 
s  =  STi/qh 

The  parameters  H,  o  ,  h,  qh,  C  ,  u^and  tu2denote  wave  height  and  fre¬ 
quency  of  incident  waves,  water  depth,  draft,  and  amplitudes  of  body 
translation  in  the  vertical  (z)  direction  and  rotations  around  x  and  y 
axes. 

The  unknown  functions  fj,  f2  ,  and  included  in  the  velocity 
potentials  satisfy  the  Helmholtz  equation  and  can  be  expressed 
respectively  by  the  Green  Functions: 

Hgl\kr),  K0(kr),  log(l/kr)  and  K0(sr) 

When  the  Green's  theorem  is  applied  to  each  sub-region  and  the 


Green  Functions  and  the  proper  boundary  conditions  at  the  artificial 
border  are  substituted  into  the  theorem,  the  equations  defining  the 
velocity  potentials  at  the  border  are  obtained.  Solving  these  equations 
along  with  the  body  motion  equations  leads  to  the  determination  of  the 
velocity  potential  at  the  border,  which,  in  turn,  can  provide  the  value 
of  velocity  potential  at  any  position  in  the  region  of  interest. 


This  model  features  the  decomposition  of  the  velocity  potential, 
which  allows  the  adoption  of  simple  form  Green  Functions  and  consequently 
results  in  economical  computation.  This  feature  also  contributes  to  the 
easy  input  data  generation  since  discretization  is  required  only  at  the 
horizontal  outline  of  the  body  because  the  vertical  distribution  func¬ 
tions,  i.e.  Eigen  Functions,  are  separated  from  the  velocity  potential. 

Garrison  Model  (GEM) 

In  this  method  the  Green  Function  is  used  which  satisfies  the 
Poisson  equation  and  all  the  boundary  conditions  but  the  one  at  the  body 
surface. 
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at  free  surface 
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The  particular  expression  for  the  Green  Function  appropriate  to  our 
boundary  value  problem  is  given  in  two  different  forms.  One  of  these  is 
the  integral  form: 
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P.  V.  denotes  the  principal  value  of  the  integral.  An  alternate  series 
form  of  the  Green  function  is  given  by: 

G  =  ~~ —  -  -  cosh  [k(C+h)]cosh  [k(z+h)][Y0  (kr )  -  iJ0(kr)] 
k2h  -  ^ht  v 


2  2 

«  Ut.  +  V 

+  — - - - cos[  uk(z+h)]  cost  Kk(C+h)  ]  K0(Uur) 

h  +  v2h  -  v  K  K  * 
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Comparing  the  Green  Functions  used  in  the  Ijima  and  Garrison  models 
it  is  readily  apparent  that  the  evaluation  of  the  Green  function  in  the 
Garrison  model  is  much  more  difficult  and  time  consuming. 

In  the  Garrison  model  the  hydrodynamic  problem  is  first  solved  and 
then  the  body  response  is  analyzed.  We  will  discuss  this  two-step  proce¬ 
dure  here.  For  simplicity  let  us  consider  the  case  of  the  one-freedcm- 
motion  system.  The  equation  of  motion  may  be  written  in  the  following 
form: 


m  +  k  a  =  0 

f  i! 

nds 

dt2 

°  -iap  i 

f*0  nd a 

(A— 9) 

where  a,  m,  k,  0,  and  So  are  the  displacement,  mass,  stiffness  due  to 
hydrostatic  and  mooring  forces,  velocity  potential  and  the  equilibrium 
immersed  body  surface,  respectively.  The  right  hand  side  of  the  equation 
represents  the  hydrodynamic  force,  which  is  only  the  drive  force  in  this 
problem.  The  velocity  potential  0  can  be  split  into  three  components  in 
this  case;  0O  associated  with  the  incident  wave  of  unit  wave  height,  0i 
associated  with  unit  body  motion  in  calm  water,  and  02  associated  with 
the  scattering  of  the  incident  wave  by  the  restrained  body. 

Because  the  velocity  potentials  due  to  incident  and  scattered  waves 
are  proportional  to  the  incident  wave  height  and  the  potential  due  to 
radiated  waves  is  linearly  related  to  the  amplitude  of  body  motion. 


A5 


the  total  velocity  potential  0  can  he  given  by: 


0  —  H  0 q  +  oi  0 1  +  II  02 


{ A-10 ) 


Substituting  this  expression  into  equation  (A-9),  we  obtain: 


m 


■^7+  ka  =  -iopH  f  (0.  +  0,  )nds  -  iapa  f  0,  nds 
dt2  JS0  */S0 

=  Fe  +  Ff  (A-ll ) 


where  Fe  is  the  force  component  due  to  the  incident  and  scattered  veloci¬ 
ty  potentials  and  Ff  is  the  component  due  to  radiated  velocity  potential 
generated  by  body  motion.  The  force  Ff,  which  is  proportional  to  the 
amplitude  of  body  motion,  may  be  conveniently  decomposed  into  two  compon¬ 
ents  in  phase  with  the  velocity  and  the  acceleration  and  may  be  written 
in  the  following  form: 


Ff 


m* 


d2u 

dt? 


(A-12 ) 


where  the  coefficients  m*  and  c*  are  termed  the  added  mass  and  damping 
coefficients,  respectively,  since  they  assume  corresponding  roles  in  the 
equation  of  motion.  Thus,  using  the  definition  of  the  hydrodynamic 
coefficients  and  rearranging  terms,  the  equation  of  motion  can  be  written 
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m  the  form: 


/  . .  ^2<ji  du 

(m  +  m*)  — j  +  c*  —  +  ka  =  F. 

at  at  c 


(A-13 ) 


In  the  above  equation,  m*,  c*  and  the  exciting  force  Fe  may  be  evaluated 
prior  to  solving  the  equation,  since  all  of  them  are  independent  of  the 
body  motion  as  discussed  above. 

Thus, in  the  first  stage, velocity  potentials  due  to  the  scattered 
waves  by  a  fixed  body  and  the  radiated  waves  by  the  unit  body  motion  in 
each  direction  are  computed,  and  these  potentials  give  the  hydrodynamic 
coefficients  m*,  c*,  and  Fe.  The  equation  of  motion,  with  all  coeffi¬ 
cients  defined,  is  then  solved. 


Appendix  B 


Subroutine  FINGFEEN 


(Extracted  Fran  "USER  MANUAL  FOR  •FINGRESN" 


By  J.  N.  NeMnsn  and  P.  D.  Sc Lavounos ) 


1 .  INTRODUCTION 

FINGREEN  is  a  file  of  subroutines  intended  for  use  in  computing  the 
velocity  potential  and  its  derivatives  for  the  3-dimensional  free  surface 
source  of  harmonic  time  dependence,  in  a  fluid  of  arbitrary  depth, 
subject  to  the  linearized  free-surface  condition.  Normally  this  file 
will  be  used  in  conjunction  with  a  radiation/diffraction  main  program 
provided  by  the  user,  which  computes  the  hydrodynamic  pressures  and 
forces  acting  upon  a  floating  or  submerged  vessel. 

FINGREEN  can  be  used  to  evaluate  either  the  infinite  or  finite-depth 
source  potential.  No  practical  limits  restrict  the  magnitude  of  the 
depth  if  it  is  finite,  nor  is  there  any  degradation  of  accuracy  if  the 
depth  is  large.  This  presents  the  user  with  a  choice  in  solving  problems 
in  relatively  deep  water,  and  comparisons  between  the  infinite-and 
finite-depth  options  may  be  made  to  determine  the  effect  of  finite  depth. 
Except  for  a  saving  in  computational  time,  there  is  no  advantage  to  using 
the  infinite-depth  option  as  an  approximation  for  large  but  finite 
depths.  The  finite-depth  option  is  recommended  in  special  cases  involv¬ 
ing  wave  interaction  with  very  deep  structures  in  close  proximity  to  the 
bottom,  where  the  fluid  is  deep  in  relation  to  the  wavelength  but  not  in 
relation  to  the  body  configuration. 

The  FINGREEN  file  consists  of  a  primary  subroutine  FGRERN  and  two 
subsidiary  subroutines  IGREFN  and  HGREFN  which  are  called  from  and  return 


*  i  ,*  >.<  • 


to  FGREKN.  This  fiLe  is  m  conventional  FORTRAN  77  source  code,  and  is 
to  be  Linked  with  the  main  program  following  compilation.  Usage  follows 
the  FORTRAN  convention,  with  a  call  from  and  return  to  the  main  program. 
Single-precision  real  numbers  are  used  throughout  the  subroutine.  The 
object  code  file  requires  17408  bytes  of  storage  on  the  VAX  11/780. 

This  manual  describes  the  use  of  FINSREEH,  provides  check  values  of 
the  output,  and  briefly  describes  the  various  parts  of  the  code. 
Reference  1  may  be  consulted  for  a  description  of  the  algorithms  which 
are  used,  and  for  definitions  of  the  source  potential  in  infinite  and 
finite  depths  (equations  1-3  and  equation  8,  respectively,  in  Reference 
1). 

2.  ACCURACY 

FINGREEN  has  been  designed  to  yield  an  accuracy  of  six  decimal 
places  for  the  source  potential  and  its  first  derivatives,  and  five 
decimal  places  for  the  second  derivitives.  Here  the  designation  of  six 
decimal  places  accuracy  is  used  to  denote  either  the  absolute  or  relative 
error,  whichever  of  these  two  is  smaller.  If  the  output  is  less  than  the 
unity  in  absolute  value  the  absolute  error  should  be  less  than  IE-06  (1E- 
05  for  the  second  derivitives)  whereas  if  the  output  is  greater  than 
unity  the  relative  error  should  be  smaller  than  the  same  amount.  This 
accuracy  test  is  applied  separately  to  each  real  component  of  the  output 
array. 

The  principal  qualification  to  this  general  statement  is  for  the 
regime  H<1,  where  the  dominant  length  scale  governing  the 
nondimensionalized  source  potential  and  its  derivatives  is  the  depth 
parameter  h,  instead  of  the  inverse  wavenumber  l/K.  For  this  reason,  the 
design  objective  in  the  regime  H<1  is  modified  to  preserve  six  decimals 
accuracy  (five  for  the  second  derivatives)  for  the  products  of  H  with 
the  potential,  H  with  the  first  derivatives,  and  H  with  the  second 


In  some  instances,  the  accuracy  may  be  degraded  by  the  accumulation 
of  round-off  errors  and/or  subtractive  cancel  lation.  One  specific  and 
obvious  example  of  subtractive  cancellation  occurs  when  the  source  and 
field  points  are  very  close  to  the  bottom  of  the  fluid  domain,  and  to 
each  other;  in  this  case  the  image  source  beneath  the  bottom  is  dominant 
(the  second  term  on  the  right  side  of  equation  10a  in  Reference  1),  and 
(since  the  origin  is  at  the  free  surface)  the  vertical  component  (2h  +  z 
+£)  involves  the  difference  between  nearly-equal  quantities. 

A  similar  problem  occurs  for  large  values  of  X,  since  the 
oscillatory  phase  of  the  output  is  limited  in  precision.  An  estimate  of 
this  error  is  the  product  of  2 n  multiplied  by  the  maximum  error  in  the 
machine  evaluations  of  sin(X*R0H)  and  cos(X*R0H).  Here  RCH  denotes  the 
nondimensional  wavenumber  corresponding  to  the  depth  H.  In  the  finite 
depth  case  this  parameter  is  limited  in  precision  to  a  minimum  of  6 
significant  figures  by  the  approximation  of  the  wavenumber  as  the  root  of 
the  dispersion  relation.  Thus  the  maximum  accuracy  of  the  outputs  is 
degraded  by  approximately  one  decimal  place  if  X=10,  or  two  decimal 
places  if  X=100.  In  the  infinite-depth  limit  (or  for  values  of  H>10, 
where  the  error  in  the  approximation  of  RDH  is  negligible)  the  limiting 
factor  is  the  accuracy  of  the  machine  evaluation  of  the  trigonometeric 
functions  of  large  argument  X,  which  is  of  similar  magnitude  to  the  above 
estimates  for  computing  systems  with  6-7  decimal  single-precision 
accuracy. 

Since  the  actual  errors  are  machine-specific,  users  who  are 
concerned  to  achieve  extreme  accuracy  are  advised  to  perform  their  own 
tests  by  comparing  the  outputs  from  FINGREEN  with  a  double-precision 
table.  Appendix  1  (Appendix  not  included)  contains  tables  of  test  out¬ 
puts  both  from  FINGREEN  (run  on  the  VAX  11/780  using  single-precision 
arithmetic)  and  from  a  special  double-precision  program  DGREFN  which  is 
designed  to  be  accurate  to  twelve  decimals.  For  the  infinite-depth  case 
double-precision  tables  of  FR  and  FR1  are  contained  in  Reference  2. 


3. 


ESTIMATES 


The  running  times  of  FINGREEN  depend  on  the  inputs,  and  the  type  of 
computing  system  which  is  used.  The  following  estimates  are  based  on  use 
of  the  VAX  11/780  system.  Limited  tests  on  an  IBM  370/168  and  an  IBM 
370/3033N  indicate  that  the  time  required  on  each  of  these  systems  is 
about  40%  and  15%  respectively  of  those  on  the  VAX.  In  other  words,  one 
call  to  FINGREEN  for  the  infinite  depth  case  which  may  require  1.0 
milliseconds  on  the  VAX  will  require  0.4  milliseconds  on  the  370/lo8,  or 
0.15  milliseconds  on  the  370/3033N.  The  same  run  on  the  IBM  PC 
microcomputer  equipped  with  a  floating  point  coprocessor  requires  about 
15  milliseconds. 

For  the  infinite-depth  case,  the  time  depends  on  the  two  coordinates 
X,  Y,  where  X=X1  has  the  same  meaning  as  in  finite  depth  and  Y  is  the  sum 
of  the  two  vertical  coordinates,  Y=X3+X6.  The  maximum  run  times  will 
occur  when  X  and  Y  are  of  comparable  magnitude;  the  maximum  time  which 
has  been  observed  on  the  VAX  is  1.2  milliseconds.  Faster  times  will 
occur  if  either  X  or  Y  is  small,  or  if  both  coordinates  are  very  large. 

For  the  finite-depth  case,  a  distincition  is  made  between  the 
initialization  of  parameters  which  depend  only  cn  H,  using  the  auxiliary 
subroutine  (GREEN,  and  the  subsequent  evaluation  in  PGREEN.  When  a  new 
value  of  H  is  input  to  FGREEN,  six  milliseconds  is  required  by  the 
initial  call  to  (GREEN.  This  step  is  avoided  for  subsequent  calls  with 
the  value  of  H  unchanged.  Thus,  in  typical  usage,  where  many  calls  are 
made  successively  with  different  sources  and  field  points  but  the  same 
value  of  H,  the  time  required  by  HGREFN  is  negligible. 

The  remaining  time  for  FGREEN  depends  primarily  on  the  ratio  X/H, 
with  a  maximum  run  time  of  6  milliseconds  occuring  when  0.5  <  X/H  <0.6. 
For  smaller  values  of  X/H,  the  run  time  will  vary  between  about  2  ms.  and 
4  ms.  as  X/H  increases  from  zero  to  0.5.  For  larger  values  of  X/H  the 
run  time  will  decrease  in  a  manner  which  is  inversely  proportional  to 


X/H.  Thus  for  many  calls  in  finite  depth  the  average  run  times  are  about 
3-4  milliseconds  cn  the  VAX. 
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